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^^ The theory of nonlinear response for Markov processes obeying a master equation is formulated in 

OO terms of time- dependent perturbation theory for the Green's functions and general expressions for 

I — I the response functions up to third order in the external field are given. The nonlinear response 

^+H is calculated for a model of dipole reorientations in an asymmetric double well potential, a stan- 

CZ3 dard model in the field of dielectric spectroscopy. The static nonlinear response is finite with the 



exception of a certain temperature Tq determined by the value of the asymmetry. In a narrow 
temperature range around Tq, the modulus of the frequency-dependent cubic response shows apeak 
at a frequency on the order of the relaxation rate and it vanishes for both, low frequencies and 



e 

C high frequencies. At temperatures at which the static response is finite (lower and higher than Tq), 

Q the modulus is found to decay monotonously from the static limit to zero at high frequencies. In 

addition, results of calculations for a trap model with a Gaussian density of states are presented. 

'""' In this case, the cubic response depends on the specific dynamical variable considered and also on 

l/^ the way the external field is coupled to the kinetics of the model. In particular, a set of different 

^ dynamical variables is considered that gives rise to identical shapes of the linear susceptibility and 

T-H only to different temperature dependencies of the relaxation times. It is found that the frequency 

cn dependence of the nonlinear response functions, however, strongly depends on the particular choice 

P^ of the variables. The results are discussed in the context of recent theoretical and experimental 

^H findings regarding the nonlinear response of supercooled liquids and glasses. 

> 

X 

I. Introduction 

In recent years progress has been achieved in the understanding of the heterogeneous dynamics 
observed in supercooled hquids and glassy systems [H [2]. Starting with NMR experiments [31 HI E] a 
number of frequency-selective techniques have been developed in order to investigate the nature of 
the dynamic heterogeneities in the slow primary relaxation of supercooled liquids[6l d El [9]. Also 
the length scale associated with the heterogeneities could be determined in some cases [TOl \TT\ . 
In the experimental studies the system always is monitored at more than two times via the 
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observation of four-time correlation functions as in the quoted NMR experiments. Alternatively, 
large external fields are applied giving rise to nonlinear effects as in the nonresonant hole burning 
studies [T2I [T3]. Furthermore, in computer simulations on model systems dynamic heterogeneities 
have been observed via following certain trajectories (THUS] or also via the calculation of four-time 
correlation functions[T6l ITT] . 

Most of the studies on dynamic heterogeneities were concerned with systems in thermal equi- 
librium, but also aging glasses have been investigated P^ ITU]. Heterogeneous aging has also been 
studied theoretically in spin glasses [2U]. in simple spin models [21] and also in a free-energy land- 
scape model for glassy relaxation [22] ■ 

In recent years, both experimental techniques and theoretical tools have been refined in order 
to allow detailed investigations of dynamic heterogeneities. In particular, it has been recognized 
that higher-order correlation functions that probe the system at different times and different 
locations in space can be used to observe a length scale [23] and the relevant four-point correlation 
function Xiit) has been studied theoretically [211 I2S1 [2S]- Earlier experimental studies used the 
approximative relation of X4(^) to a two-point correlation function [23| |27] in order to extract the 
number of cooperatively rearranging particles, A^corr- In an influential paper Bouchaud and Biroli 
related the nonlinear (cubic) response Xai^^T) to X4(^)[28]. The experimental determination of 
Xsi^JjT) allowed the determination A^corr more directly [211130] and the results are compatible with 
the earlier observations. In particular, it was argued that the function 

XM = \x:MT)\^, (1) 

with Axi denoting the static linear response, ks the Boltzmann constant and a^ the molecular 
volume, exhibits a hump-like structure. This behavior is assumed to be a distinctive feature 
of glassy correlations [21]. Additionally, the maximum of X{uj,T) is expected to decrease with 
increasing temperature and to be directly proportional to Acorr- If glassy correlations are absent, 
X{u},T) should not be peaked and this 'trivial' behavior consists in a smooth cross-over from 
a low-frequency limiting value to a vanishing high-frequency limit. In this context it has to be 
mentioned that Brun et al. found a hump-like shape for X{u,T) in a calculation employing the 
so-called box model[3T], a model devoid of spatial aspects. 

Apart from the determination of A"corr. the nonlinear dielectric response has been used to 
investigate the nature of the heterogenous dynamics via comparison of the cubic response with 
the linear response [321 EH] and the results were discussed in the framework of the box model. 
Similar measurements were performed in order to extract the configurational heat capacity of 
liquids[3l|. In addition also the nonlinear dielectric response of liquids due to an AC and an DC 
field pulse have been recorded[35| and also dipolar glasses have been investigated [36]. 

The present paper deals with the theory of nonlinear response functions for Markov processes, 
because the relaxation in complex systems often is modeled in terms of such stochastic dynamics. 
For systems that follow a Hamiltonian or Langevin dynamics, nonlinear response functions have 
been considered quite some time ago [371 l38l 139] . However, explicit calculations of response func- 
tions are rare and most of them relate to variants of the rotational diffusion of molecules in the 
presence of strong electric fields, see e.g. refs. [TO FiTl H2 l H3] . In addition, approximate nonlinear 
response theory has been investigated more generally [ii] and also fluctuation-dissipation relations 



beyond the linear regime have been discussed [l5l H6] . The nonhnear response of supercooled liq- 
uids has been worked out theoretically in the framework of mode-coupling theory |17j. Here, I will 
perform the calculation of the response functions in close analogy to the quantum-mechanical way 
of computing response functions [48j . Time-dependent perturbation theory for the propagator is 
used in order to obtain the response in the desired order in the amplitude of the external field. I 
will present the results of calculations of the cubic response function for two Markovian models 
of relaxation. One model describes the reorientations of dipoles in an asymmetric double well 
potential (ADWP) and has been used to interpret results of dielectric experiments in general [lU]. 
Furthermore, it has also been employed in calculations of the signals obtained in nonresonant hole- 
buring experiments [50]. It will be shown that X{u,T) mainly behaves 'trivially' for this model. 
Another model that will be considered is the trap model with a Gaussian density of states [HIl E2]- 
This model has been used in the interpretation of some features of the relaxation in simulated 
supercooled liquids, both in equilibrium [33] and in the aging regime [SH ES]. Here, the results for 
X{u,T) are more complex and, depending on the parameters chosen, either exhibit a peak-like 
structure or 'trivial' behavior. 

The paper is organized as follows. In the next section, I will outline the calculation of nonlinear 
response functions for systems obeying a master equation. For convenience of the reader, most 
of the explicit calculations are presented in the Appendix. The sections following this theoretical 
part deal with a discussion of the results obtained for the two models considered and the paper 
closes with some concluding remarks. 

II. Nonlinear response theory for Markov processes 

In this section, I will outline the general procedure to calculate the nonlinear response functions 
for a system that is described by a master equation (ME) [561 EZ] • If one is dealing with complex 
systems a coarse-grained procedure may result in a description of the underlying dynamics in 
terms of a non-stationary Markov process. Therefore, in order to keep the treatment general, I 
will treat the case of a ME with time-dependent transition rates. 

In the following, Gki{t,to) denotes the conditional probability to find the system in state k at 
time t provided it was in state / at time to (Green's function, propagator) in a discrete notation. 
If continuous variables are considered, all sums in the following expressions are to be replaced by 
the corresponding integrals. Denoting the rates for a transition from state k to state I by Wik{t), 
the ME reads: 

^Gklit, to) = - E Wnk{t)Gkl{t, to) + J2 Wkn{t)Gnl{t, to) (2) 

This equation has to be solved with the initial condition Gkiito^to) = 6ki, where 6ki denotes the 
Kronecker symbol. If the transition rates Wki{t) are time-independent the process considered is 
stationary. The one-time probabilities Pk{t) (the populations of the states) obey the same ME 
and are given by pk(t) = J2iGki{t,tQ)pi{to). The Wkiit) can be related to the elements of the 
master-operator yV{t) via|56]: 

W(t),i = Wki{t) - 6ki E Wniit) (3) 

n 



Here, W(t)fc/ > holds for all A; 7^ / and the sum rule X!fc VV(t)fc« = is fulfilled for all values of 
/ as it is a general property of the transition rates for any Markov process. At the initial time to 
the system is described by a fixed set of populations, p° =Pk(to) with Y^kPk = 1- If ^ stationary 
system is considered, one often starts from equilibrium populations p^ = p^*^ or if one is interested 
in describing a situation with a certain thermal history one might choose the p^ as the equilibrium 
populations at a temperature different from the working temperature. 

In order to treat the system in the presence of an external field one has to specify the field- 
dependence of the transition rates, which is not straightforward. In case of Hamiltonian or 
Langevin dynamics, the linear coupling of a variable M{t) to a field H{t) gives rise to an ex- 
tra term [— M (t) ■ H{t)] in the Hamiltonian. In a Fokker-Planck equation, this gives rise to a term 
linear in if jSSj. If one considers a ME, one choice that has been used in a number of investigations 
of fiuctuation-disspiation relations is given by 

Wif\t) = W^H(t)e'^^[^^^^-^"^'l (4) 

with arbitrary 7 and /i[59l IMl El]. In this expression (3 = T~^ denotes the inverse temperature 
with the Boltzmann constant set to unity, ks = 1- If the system obeys detailed balance, one has 
the restriction 7 + /i = 1. In particular, for systems described by a Fokker-Planck equation, one 
would naturally choose 7 = /i = 1/2 and a linear expansion of eq.Q gives the usual term in the 
Fokker-Planck operator. However, it is obvious from eq.Q that in general one will have nonlinear 
contributions to the perturbation also if the coupling to the field is linear in the sense described 
above. This means that couplings of a form like [M{t) ■ H"^], as it would appear for instance if the 
couphng to an induced dipole-moment is considered PU]. are absent. 

In order to keep the treatment general, I will formulate the response theory without fixing the 
field-dependence of the transition rates. It is only assumed that it can be cast in the form: 



00 1 Jn 

wif\t) = j:^wj^\t).[pH{t)r with wi^\t) = ^^^^^wtf\t) 



(5) 

The elements of the propagator G' (t, to) are obtained from the ME, eq.O), where the field- 
independent quantities are replaced by those explicitly depending on the external field, i.e. 
Gif (t,to) = -En W^if (t)4f (t,to) + j:nWl:^\t)Gif{t,to). The solution of this equation is 
needed to calculate the response of the system to an external field applied at time to and mea- 
sured by an observable F(t), 

(^W)w=E^^4f(^'^oK(to) (6) 

kl 

In order to be able to set up a perturbation theory for G'^^^(t,to) in terms of the corresponding 
'field-free' propagator G(t,to), one uses the decomposition 



W(^)(t) = >V(t) + V(t) with V(t) = ^V(")(^) (7) 

where the perturbation is given according to eq.([5]) 



n=l 



n! 



wi?\t)-Skij:wif{t) 



(8) 



The theoretical treatment is very similar to the one utilized in ref.[6T] and consists in performing 
time-dependent perturbation theory to treat V(t) in the desired order of the field. The details of 
this procedure are described in Appendix A. The explicit expressions for the response functions 
are given up to third order in the field and the extention to higher order is straightforward. 

The main difference to the formalism utilized for Hamiltonian or Langevin dynamics with a 
linear coupling to the external field is that here in general the elements V^'^\t)ki with n > 1 do 
not vanish. This gives rise to a number of extra terms. The situation is visualized in FigjT| which 
shows the diagrams representing the interaction with the field for the third-order response. One 
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Figure 1: Pictorial representation of the perturbation expansion for the third-order response. The 
unperturbed propagators are denoted by G and the V*^"-* are the perturbations according to eq.(p|. 



has the terms stemming from purely linear interactions given in the first line. These terms also 
appear in a Fokker-Planck treatment of a linear coupling. Furthermore, one has two cross terms 
between first-order and second-order perturbations (second and third line in FigjI]) and a term 
stemming from the third-order perturbation (fourth line). For Langevin dynamics, cross-terms 
only appear if a quadratic coupling is considered in addition to a linear one. 

While in Appendix A the general expressions for the response functions are given, in the actual 
model calculations I will consider only the response of systems that are in thermal equilibrium 
prior to the application of the external field. Furthermore, the models treated in the present paper 
represent stationary Markov processes with time-independent transition rates. The discussion will 
be limited to sinusoidal fields of the form 

H{t) = Ho cos (oot) (9) 

For this oscillating field the linear and the cubic response for times long compared to the initial 
transients can be written as: 



Ho 
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^W(t) = _^ e— *xi(^) + c.c. 



-iuit^(l){, ,\ I ^-i3ujt^(3)i 



^(3)(i) = _o e-^^'x's'iu) + e-^'^'x's'iu;) + cc. (10) 



where c.c. denotes the complex conjugate. 

In the following sections, I will mainly discuss the quantity X(ci;,T) introduced in eq.Q. As 
the models that will be considered in the following are not related to any spatial aspects of dipole 
reorientations or relaxing units, the molecular volume will be set to unity, a^ = 1. Additionally, 
one has a separate function for each frequency-component, cf. ref. [30], that can be written as 
(a = 1,3): 

X^{uj,T) = ^^^\xt\uJ,T)\ (11) 

Thi s fun ction eliminates the 'trivial' te mper ature dependence of Xz V^i"^) because Axi ~ /3 , cf. 
eq.(A.6) and Xs ~ Z?'^ according to eq.(A.9). Therefore, any temperature dependence stems from 



the 'intrinsic' relaxation behavior of the dynamical variable considered. 

III. The ADWP-model for dipole reorientations 

In this section, I will present the results for one of the simplest models for dielectric relaxation, 
namely the model of dipole reorientation in an asymmetric double well potential. I will closely 
follow the notation used in a related investigation of the nonresonant dielectric hole burning 

technique da [131 EH]. 

As in ref.jnO], two dipole orientations denoted by '1' and '2', characterized by polar angles 
9i = 9 and 6*2 = 6' + vr are assumed and the transition rates between the two are given by 
Wi2 = We~^^^'^ and W21 = We'^^^/'^. Here A denotes the asymmetry, and W is the hopping rate 
in the symmetric case. For this model, the Green's functions in the field-free case are are given 
by: 

Gkiit) = pf (1 - e"*/") + 5kie-''^ with r"^ = 2W cosh(/3A/2) and p"^ = t ■ Wm (12) 

The variable that couples to the field is 

Mk = M cos{9k) and therefore Mi = M cos{9) ; M2 = -M cos(6') 

with M denoting the static molecular dipole moment. The field- dependent transition rates are 
chosen as in eq.(|4]) with 7 = /i = 1/2. (If this restriction is relaxed all response functions depend 
on the sum (7 + yu), which equals unity in the present case.) In the calculation of the response 
I assume a collection of systems characterized by an isotropic distribution of orientations and 
therefore an average over the angle 9 is performed according to (cos"(^fc)) = (n + 1)~^ for n even 
and (cos"(6'a;)) = for n odd. 



Using the general expressions given in Appendix A along with eq.(12), one finds for the linear 
response: 

XiM = Axi^^— where Axi = P{AM') = ^^ (l ~ 6') (13) 

1 — tUJT 3 ^ -^ 

In this expression, I defined 6 = tanh(/3A/2). (It should be mentioned that Axi differs by 
a factor 1/2 from the definition of Xdwp in ref. |50^.) As usual, Axi is related to the mean- 
square fluctuations of the dipole moment (AM^). Eq.(13) follows immediately from the definition 



{M"^) = Y^k^k^P^k ^^d eq.(12) with additional isotropic average. 
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Note that in the ADWP-model, the static susceptibihty A^i for non-vanishing asymmetry de- 
pends on temperature due to the dependence on 6 in addition to the trivial 1/T-dependence. This 
behavior for finite asymmetry is different from the model of Brownian rotational diffusion, where 
TAxi is independent of temperature [19]. For vanishing asymmetry, the models show identical 
behavior (apart from irrelevant prefactors). Without showing results here, it is mentioned that 
Re(xi(w)) decays from its low-frequency limit A^i to zero for large frequencies and Im(xi(a;)) 
shows the typical Lorentzian behavior and is peaked at ut = 1. 



The third-order response functions are calculated according to eq.(lO) using the general ex- 



pressions given in eq.(A.9) in the Appendix. In a straightforward calculation one finds: 



M 



X^M = 9-/3^(1- 5^) x5r^(c.r) 



20 
Here, the spectral functions only depend on the product x = ut and are given by: 



:(i) 



S's>{x) = 6'' 



3{l + i2x) 



+ 



2{x^ - 1) + ix{x^ - 3) 



2(1 



X 



2^2 



Si'\x) 



;i + x2)(i + 4x2) 

(1 - 11x2) + i6x{l - x^) 2(5a;2 - 1) + i3x{x^ - 3) 



(14) 



(15) 



;i + x2)(l + 4x2)(l + 9a;2; 



+ 



6(l + x2)(l + 9x2) 



When compared to the model of Brownian rotational diffusion, the following can be observed. For 



.i^)i 






A = 0, Xs {^) ioT the two models are very similar, cf. FigJ2 and Figs. 3,4 of ref. j41j. For finite 
A, however, the third-order response for the ADWP-model shows a characteristic temperature 
dependence, that is absent in the model of rotational Brownian motion. 

In FigM the real and the imaginary part of the 30;- component Xs {^) ^ire plotted versus cot 
for different values of the asymmetry A and various temperatures. It is evident that the sign of 
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Figure 2: Real part (red) and imaginary part (black) of the 3a;-component Xs (i^) for the ADWP- 



model as a function of wr, where r is the relaxation time according to eq.(12). 

both functions change as a function of frequency. Furthermore, the shapes of Im(x3 (cu)) differ 



7 



significantly from Lorentzians. As mentioned above, for A = 0, Xs i^) does not depend on 
temperature. 

Tlie static nonlinear susceptibilites are determined by the limiting values of the spectral dfunc- 
tions, Si^\0) = (352 - 1) and Si^\o) = {35'^ - l)/3, and thus are given by: 



.(3) 



x's'iO) 



1)0 



13'' (352 



1-5' 



.(1) 



xr(o) = 3xr(o) 



(3), 



(16) 



It should be mentioned, that X3 (0) is determined by the fourth-order cumulant, K,i{M) = (M^) — 
4(M)(M3) - 3(M2)2 + 12(M)2(M2) - 6(M)^ = 2M^ {36^ - 1) (1 - 5^)- For finite A, the low- 
frequency limit X3 (0) vanishes at a temperature Tq, at which 5*3 (0) = 0, 

To = A/ In [(v^ + l)/( VS - 1)] ~ A/1.317. 

For large frequencies, one always has X3 (00) = 0. 

Instead of discussing Xs i^) further, in the following I will consider Xa{uj,T) according to 



eq.(ll). This quantity is given by, cf. eq.(13) and eq.(|14|): 
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XaicO,T) 



:(") 



sn^r) 



20 (1-52) 



(17) 



The limiting values for small and large frequencies are determined by the corresponding limits of 
5*3 (wr) and thus, one has for example X3(0,T) = {3/20){\36'^ — 1| / (1 — 5'^)). It is evident, that 
Xa{uj, T) will have a peak-like structure for T ^ Tq. As is shown in FigM for other temperatures 
one has 'trivial' behavior, i.e. a continuous decay from the low-frequency limit to Xq(cj,T) = 
at high frequencies. One can see, that the behavior of the Iw-component and the 3a;-component 
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Figure 3: Xq,(ci;, T) for various values of the asymmetry and different temperatures. In the upper- 



most panel, X 



(Debye) , 



a;)[lT] is shown for comparison (dashed line). 



is very similar. 
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Figure 4: X™^^(u;) 7X3(0) versus temperatures for A = 1. The dotted line is the same with 
assumption of a Gaussian distribution of A with mean A = 1 and variance cta = 10. 

In order to further quantify the behavior of Xa{u},T) with regard to a 'hump'-hke structure, 
in Fig|4| the ratio X^^^{u)/X^{0) is plotted versus temperature. For T ^ Tq and also for T ^ Tq 
trivial behavior is observed and only in the region of T ~ Tq a hump develops. This hump, 
however, has nothing to do with glassy correlations but is solely a consequence of the temperature 
dependence of the fluctuations of the dipole moments. 

Finally, it is to be mentioned that the above results hardly change if one considers distributions 
of the hopping rate W and/or the asymmetry. In particular, the temperature-dependent change 
in the shape of Xa{co) is practically unaltered. This is exemplified in Fig El where the dotted line 
represents X™'*^(w) 7X3(0) for the case of a broad Gaussian distribution of A. The reason for this 
is simply the steepness of the root of 6*3 (0) = 0, meaning that the overall behavior is determined 
by the mean value of A. Thus, if one considers a system with a distribution of asymmetries that 
is centered at A = 0, one will observe trivial behavior of Xa{uj) at all temperatures. Ladieu et al. 
use the ADWP-model with finite A and some further assumptions to fit the experimental data on 
supercooled liquids [62]. 



IV. Trap models 

In this section, I will discuss Xa{co, T) for the trap model with a Gaussian density of states, which, 
as already mentioned in the Introduction, shows some features of glassy relaxation. It is defined 
by the ME for G(e, t + to|eo, to) = G{e, t\eo, 0) = G{e, t|eo), in a continuous form written as: 

^(e, t|eo) = -K(e)G'(e, t|eo) + p(e) ( de' K{e')G{t\t\eo) (18) 



In eq.(18), the escape rate is given by 



Kiel — KinrtG' 



P^ 



(19) 



with the attempt rate n^o. Furthermore, I solely consider the model with a Gaussian DOS 



p{e) 



J_^-.2/(2<x2) 



27rc7 



(20) 



with cr = 1. From eq.(18), the equihbriuin populations at a given temperature T (measured in 
units of 0") are found to be Gaussian p^'^(e) = hmj_j.oo G{e, t\eo) = ~7^^~^''~''^ ''^^'^ ^ with e = —(3a'^. 
In order to calculate the response, one further has to quantify the dependence of the dynamical 
variable on the trap energy e. The choice of this dependence represents a further assumption of the 
calculation and has a strong impact on the results for the cubic response, as will be discussed below. 
In order to clarify this issue, consider the linear response for the specific choice of eq.(|4]) for the 



field-dependence of the transition rates. Using eqns.([6]), (A. 5) and (A. 6), one obtains the relation 



between the linear response and the equilibrium auto-correlation function CM{t) = {M{t)M{0)), 
RjJ{t) = — /3(7 + ^)[dCM{t)/dt], if the system is in thermal equilibrium [ni]. In the frequency- 
domain, this yields eq.(|B.2) in Appendix B, if the average over the possible realizations of the 



variables is performed with the following assumption: 

(M(e)) = and {M{e)M{e^)) = 6{t - t^){M{tf) (21) 

In the calculation of the third-order response, the fourth moments of the variable are important. 
For the corresponding averages I will assume a Gaussian factorization property for simplicity: 

(M(ei)M(e2)M(e3)M(e4)) = 5(ei - t2)6{t^ - t^){M{eif){M{t:,f) 

+ 5(ei - 63)^(62 - e4)(M(ei)2)(M(e2)2) (22) 

+ 5(ei - 64)^(62 - e3)(M(ei)2)(M(e2)') 

In the calculation of the response, the field-dependence of the transition rates has to be fixed 
additionally. I use eq.(|4]) with arbitrary values for 7 and /i. From the physics of the model one 
might argue that /i = 1 and 7 = is an appropriate choice because it is meaningful to assume that 
the activation energy of the escape is biased by the field (according to e — )■ e — M(e) ■ H). However, 
it is not clear that this simple argument holds in out-of-equilibrium situations and for strong fields. 
Using the assumptions made, one can compute the response according to the expressions given in 
Appendix A. The calculation is outlined in Appendix B and here only the results will be discussed. 
In the explicit choice of the variable, I follow Fielding and Sollich[63] and use a set of variables 
with an Arrhenius-like dependence on the trap energies: 

(M(e)2) = e-"''^ (23) 

with variable n and where the static value of M^ has been set to unity. For n = 0, one has 
temperature-independent variables as in case of Brownian rotational diffusion. 



The most important consequence of the specific choice, eq.(23), is that it does not affect the 
spectral shape of the linear response. The only quantities that strongly depend on the choice of 
n are the static susceptibility and the the temperature dependence of the relaxation time. This is 
because one can write: 



M'"-*^(^-^^'"M: 



eg ^lej 



K{e) - iujn 
with 



q2„2 



ujn = uje"" " (24) 
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Thus, the susceptibihty is given by: 

X^ioj) = /3(7 + /i) fdep{ere~-^^ f^^ = Axi f dep{er . ".^'^ (25) 

J K[e) — lUJ J K[e) — tUJn 

The static susceptibhty, i.e. the amphtude, Axi, strongly depends on the choice of n and reads 



as: 



Axi = (7 + f^)P{M')T = (7 + /x)/3e^'^^''^' 



(26) 



Here, the second moment {M'^)rp is related to the low- frequency limit of xi{^)^ (^^)r ~ 
Jde{M{ef')p{eY'^. Note that Axi is temperature independent only for n = and for n = —2. 
In Figjsl the imaginary part of Xi{^) is shown for n = and various temperatures. The 




10"' 10"' 10° 10' 10' 10^ 



COT. 



eq 



Figure 5: Imaginary part of Txi{(jj), Tx"(a;), for n = and various temperatures (T/cr=0.5, 0.6, 
0.7, 0.8, 0.9, 1.0 as indicated by the arrow). The dotted line represents a Lorentzian. 



frequencies are scaled to the relaxation time of Cuif) for n = 0, Teq = j^dtCuif) = '^^^^'^ '^ , cf. 
ref. [HI]. It is obvious that Xi('^) broadens as temperature is decreased and thus time-temperature- 
superposition is not obeyed. It is stressed again, that xi {^) is basically independent of the choice 
of n. 

Next, the behavior of the cubic response and its dependence on the model parameters will be 
discussed. Using the limiting values of the cubic response functions given in Appendix B for small 
and high frequencies, one finds the following limits for Xs i^)'- 

1 



.(3) 



Xr(0) = -/3^(7 + /i)'(6-6) ; X^^^(0) = 3x^(0) and xr{oo) = 



Wi 



M)l 



S»)i 



(27) 



Here, I defined the averages ^i = {M'^) ^{M"^) j, and ,^2 = {M'^^rp, which for the Gaussian trap 
model are given by: 

^^ = e«("+l)/3V2 . ^^ ^ g2n(n+l)/32a2 ^28) 

With these quantities, one finds for the low-frequency limit of X^: 

X3(0,T) = ^(7 + /i)-'^'~^'' 



{{M%y 



(29) 



and similarly for Xi(0, T). It is thus clear that these low-frequency limits do strongly depend on 
the variable, i.e. on n. Therefore, one can expect to find trivial or hump-like behavior of Xa{uj, T), 
a = 1,3, depending on this choice. 
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In FigJ6^X3(a;, T) is shown for n = and various values of /i. Here, it is assumed that 7+/i = 1. 
The main difference between the various choices for fi is the overall amplitude. Additionally, it 
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eq 

Figure 6: X3(ci;,T) for n = and various values of n for 7 = 1 — yu and different temperatures 
(T/cr = 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1) in the order indicated by the arrow. 



is clear that X^lujT) exhibits a hump in all cases. However, in contrast to the results obtained 
on supercooled liquids, the maximum value of X3 increases as a function of temperature. This 
increase is somewhat stronger for n = 1 than it is for other values of /x. 

Next, I will consider values for n different from zero, meaning that the dynamical variable that 
couples to field shows an explicit dependence on the trap energies. In FigJT^, X^(uj) is plotted 
versus frequency for n = 1 and the same values for /i as in Fig|6} It is observed that a hump 
is found at high temperatures, whereas trivial behavior is observed at low temperatures. The 
temperature, at which a visible peak is observed depends on the value of /i, i.e. on the way, 
the field couples to the transition rates. This is shown in FigjTJD, where the maximum value of 
X^lij) is plotted versus temperature for temperatures higher than the onset temperature, which 
is defined by the first appearance of a hump in Xs(u) indicated by the dots in FigJTJD. In the 
temperature range of a hump-like shape of X^{u,T) its maximum, X^'^^{u,T), appears to be 
almost independent of temperature. A similar behavior is found for other positive values of n. 

From these model calculations it becomes apparent that the existence of a hump depends on 
the value of Xa{0), the value of the maximum of Xa{uj), and in particular their ratio. Thus, the 
low- frequency limit plays an important role in determining the overall shape of Xa{u}). 

These considerations can be further substantiated by considering the special value of ra = —1, 
because in this case one has ^1 = ^2 = and therefore Xa{0) = 0, cf. eq.(28). This means, a 
hump will be observed in this case, as is confirmed in Figjsk, where Xs^co) is plotted as a function 
of frequency for /i = 1. For other values of /i, the results are very similar. On first sight, the 
behavior of Xs^co) is very similar to that for n = 0, cf. FigJGJ However, the maximum for n = —1, 



X^^^{u, T), is a decreasing function of temperature as opposed to the case of n = 0, cf. FigjSJD. 
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Figure 7: a: (left) X3{u,T) for n = 1 and different temperatures (T/a = 1,1.5,2,2.5,3). The 
arrow indicates increasing temperature, b: (right) X^^^{uj,T) as a function of temperature for 
n = 1. The curves are shown for temperatures higher than the onset temperature, below which 

0. The dotted line is the result for n = 0, /x = 1. 



trivial behavior is observed, i.e. Wmax 
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Figure 8: a: (left) X3{u, T) for n = —1, /i = 1 and different temperatures (T/a = 0.6, 0.7, 0.8, 0.9, 1 
from top to bottom), b: (right) X^'^^{u,T) as a function of temperature for n = —1 and n = 
(/x = l). 



At this point, however, it has to be noted that the case ra = — 1 is somewhat special as the mean 
relaxation time of the linear response, (r) = J dep{eY'^{e~'^'' / K{e)) = k^, is basically temperature- 
independent. Thus, although the shapes of Xi{^) ^-^e identical for n = and n = —1 at a given 
temperature, the mean relaxation time for n = —1 does not change with temperature. This shows 
that it is not straightforward to compare linear and nonlinear response functions. 
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V. Conclusions 

The theory of nonhnear response functions for a system obeying a master equation has been 
formulated in close analogy to quantum mechanical nonlinear response theory. Time-dependent 
perturbation theory is used in order to compute the elements of the propagator (the Green's 
function or conditional probability) in the desired order of the amplitude of the applied external 
field. Expressions for the response functions up to third order are given in terms of the solution 
of the field-free master equation for systems with arbitrary initial conditions and also for non- 
stationary Markov processes. In the actual model calculations, however, only stationary systems 
are considered that were in thermal equilibrium prior to the application of the field. The treatment 
of aging systems or other non-equilibrium situations are beyond the scope of the present paper. 

For the model of dipole reorientations in an asymmetric double well potential (ADWP-model), 
the spectral shape of the modulus of the frequency- dependent cubic response, X3(co',T), shows a 
specific temperature dependence which strongly depends on the value of the static susceptibility, 
^3(0, T). At a temperature Tq, which is determined by the value of the asymmetry of the poten- 
tial, ^3(0, Tq) vanishes. For a narrow temperature range in the vicinity of Tq a peak is observed in 
the modulus. For temperatures sufficiently different from Tq a monotonous decay from X3(0) 7^ 
to X3(oo) = is found. This 'trivial' behavior is basically the same as for the model of rota- 
tional Brownian motion [H] and is at variance with experimental results obtained for supercooled 
glycerol[2niEn]- It was attributed to trivial dipole reorientations that occur independent of glassy 
correlations. These correlations should give rise to a peaked behavior, i.e. the existence of a hump 
in X3(co'). If one intends to utilize the ADWP model for the dipole reorientations in supercooled 
liquids, it is natural to assume distributions of relaxation rates and of asymmetries. However, as 
shown in Section III, such a distribution hardly affects the spectral shape of X^lu) (apart from 
the fact that a distribution of relaxation times gives rise to a broadening). 

If a trap model with a Gaussian distribution of trap energies is considered, a more complex 
dependence of X^^u) on the parameters used in the calculations is observed. In particular, the 
dependence of the dynamical variables that couple to the external field on the trap energies, 
M(e), has to be fixed. I restricted the calculations to variables that obey Gaussian statistics 
and depend on the trap energy in an exponential way, M(e) = e""^"^, cf. eq.(23). This choice 



is particularly useful when discussing the properties of nonlinear response functions and their 
relation to the linear response. This is because the exponential dependence on the trap energies 
has the interesting property that the spectral shape of the linear susceptibility is the same for all 
values of the parameter n. Only the amplitude (Axi) and the temperature dependence of the 
relaxation time strongly depend on its specific value. If the nonlinear response is considered, it 
is however found that the temperature-dependent spectral shape of X^^u) strongly depends on 
the value of n. In particular, one can find a peak or 'trivial' behavior depending on both, the 
value of n and the temperature. Similar to the situation in the ADWP-model, the existence of 
a peak is related to the value of the static susceptibility. In case of the occurrence of a hump, 
the temperature dependence of the peak maximum, X^'^^lu), can increase (n = 0) or decrease 
(n = —1) with increasing temperature. These results indicate that it is difficult to compare the 
linear and nonlinear susceptibilities. It is left for future work to investigate the behavior of the 
nonlinear response in the trap model for other dynamical variables and also for non-equilibrium 
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situations. 

In the experimental determination of X3(a;,T) in supercooled liquids[29l [30], the decrease of 
Xf^'^^{u) with increasing temperature has been used to extract the number of correlated molecules, 
-^corr, which is a 'real space property' of the dynamical heterogenities in glasses. Due to the 
mean-field nature of both models considered in the present paper, none of the results presented 
have any connection to real space. Therefore, a direct comparison to experimental data is not 
possible. However, the model calculations substantiate the fact observed earlier already pi] that 
the existence of a peak in X^lu) does not have to be related to glassy correlations in some sense. 

In conclusion, I have formulated a theory of nonlinear response for systems described by Markov 
processes and have presented the results of calculations for simple stochastic models. The most 
important result is that the spectral shape of the nonlinear (cubic) response can vary considerably 
depending on the model considered. The occurrence of a peak in the modulus of the third-order 
susceptibility cannot generally be attributed to glassy correlations. Of course, this does not mean 
that glassy correlations do not give rise to a hump but its mere existence cannot be taken as 
a signature of such correlations. Because the models considered in the present paper are of a 
mean-field nature, it is impossible to connect the results to a length scale of any kind. Due to the 
growing interest in nonlinear responses in complex systems, calculations of the kind presented in 
the present paper should be performed for a variety of different models in order to gain a deeper 
understanding of the general features governing the shape and the temperature dependence of the 
corresponding susceptibilities. 

Acknowledgment 

I thank Roland Bohmer, Gerald Hinze, Francois Ladieu and Jeppe Dyre for fruitful discussions 
and Roland Bohmer for helpful comments on the mansucript. 



15 



Appendix A: Calculation of nonlinear response functions 

In this appendix the calculation of the response for a system obeying the ME, eq.([2]), using time- 
dependent perturbation theory is described. Using eq.([3]) for the master-operator, the ME in a 
matrix notation reads: 

9iG(t,to) = W(t)G(t,to) (A.l) 

Here, the propagator has matrix-elements G(t,to)w = Gki{t,to). The solution of the ME in the 
absence of an external field can be written in the form: 

G(t,to) = Texp I^J' dTW{T)^G{to,to) (A.2) 

where T denotes the time-ordering operator and G{tQ,to)ki = Ski- In the presence of the field the 
transition rates are given by eq.([5]) and the corresponding master-operator accordingly reads as 
>V(^)(t)fcz = wiPit) - dMEnW^it). The ME is written as 9tG(^)(t,to) = W(^)(t)G(^)(t,to). 

In order to calculate the response of the system to an external field applied at time t = 
to and measured by an observable F(t), {F(t))(_H) = Y^kiFkGki (t,to)pk(to) as given in eq.(6j), 
time-dependent perturbation theory is used to express the propagator as a series of the form 
G(-f^)(t,to) = G(t,to) + E^=iG(")(t,to), where G(t,to) denotes the propagator in the field-free 
case. In order to perform the calculation, one proceeds in the following way. Starting from the 
Dyson-like equation 



G^''\t,to) = G(t,to) + / dt'G{t,t')V{t')G'^''\t',to) (A.3) 

Jto 

one obtains, using eq.Q, for the lowest order terms: 
G^'\t,to)= fdt'G{t,t')V^'\t')G{t',to) 

Jto 

G(2)(t,to) = f'dt'G{t,t')V^^\t')G{t',to) + f'dt'G{t,t')V'-^\t')G^^\t',to) 

Jto Jto 

G(3)(t,to)= f'dt'G{t,t')V'-^\t')G{t',to)+ t dt'G{t,t')V^^\t')G^^\t\U) (A.4) 

Jto Jto 

+ tdt'G{t,t')V^^\t')G^^\t\to) 

Jto 

In the next step, one uses the expression for the matrix elements of G^"''(t,to); denoted by 
Gj^i {t,to), in eq.(|6) in order to compute the nth-order response, Xf (''^!^o)• 
With the definition 

4f (^2,tl) = E [Gkm{t2,h) - Gk,{t2,h)] W!^]{h) (A.5) 

m 

where Wj^j{ti) is given in eq.(5), one obtains in a straightforward calculation for the linear re- 
sponse: 

x'F\t,to)= fdhH{h)R?{t,h) with R^j^\t,h) = PY.FkL'i}i>{t,h)pi{h) (A.6) 

•^*o k,l 
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From the structure of this expression it is evident that R denotes the usual response to a short 
field kick. 

The second-order response is found to consist of two terms: 



xP{t,to) = xf'\t,to) + xP'\t,to) (A.7) 



with 



Xp'\t,to)= f dhH{h) r dt2H{t2)Rf'^\tMM) 

Jto Jto 



k,l,m 

Xp'\t,to) = ^ fdt,H{t,fRp'\tM) (A.8) 

Rp'\t,t,) = P'Y.FkL'iht,h)piih) 

k,l 

This second order response is expected to be of little relevance in most cases as it vanishes in 
isotropic systems. More interesting is the third-order response because usually this is the lowest- 
order nonlinear contribution to the response of the system. As can be expected from FigjT| it has 
the form: 

xPit.to) = X?''\t,to) + xf''\t,to) + xf'\t,to) (A.9) 

and the individual terms are given by: 

xf'\t,to) = fdtiH{ti) rdt2H{t2) r dhH{U)Rf^\tMMM) (A.IO) 

Jto Jto Jto 

k,l,m,n 

originating from the linear perturbation. This term also is found in the response theory for a 
Fokker-Planck equation. The cross-terms between the first- and second-order perturbations are: 

Xp'\t,to) = Xp'^\t,to) + Xp'''\t,to) 

xT^\tM) = \ fdhH{ti) rdt2H{t2fRp^^\tMM) 
Z Jto Jto 

Rp^^\tMM) = P'Y: FkL^klit,h)L'^^l{tuh)Mt2) (A.ll) 

k,l,m 

Xp'''\t,to) = I fdUH{uf rdhHih)Rf'''\t,h,h) 

Z Jto Jto 

k,l,m 

Finally, the third-order contribution is: 

Xp'\t,to) = ^ fdUH{ufRf'\tM) (A.12) 

Jto 

k,l 
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These expressions are valid for any Markov process obeying the ME, eq.([2]) and arbitrary initial 
conditions (initial populations p; (to)). 

If the process considered is stationary, meaning that the transition rates are time- independent, 
Wkiit) = Wkh the Green's functions depend only on the time-differences, Gki{t2, h) = Gki{t2 — h). 
Furthermore, if the system was in equilibrium initially, pi{to) = p^^, the expressions simplify 
considerably. In this case, the integrals can easily be transformed in order to find expressions for 
x''"^(t — to) that are reminiscent of the standard ones used for instance in the field of nonlinear 
opticsjlS]. If one is interested in the stationary response, one just starts recording x^^K^ — io) 
after times long compared to the initial transients. 



Finally, the explicit choice of the field-dependence of the transition rates enters via eq.(A.5). 



Appendix B: Nonlinear response functions for the trap 
model 

Using the general expressions given in Appendix A, one can calculate the response for the trap 
model. In contrast to other models, the calculation is simplified by the fact that for large N, the 
number of states, one has to consider only terms of order unity and one can neglect all terms 
of order 1/A^. In a discrete notation, one has for the trap model Gki{t) = (5fc;e~'*'=* + 0{1/N) 
with Gkiif) = G{ek,t\ei) and n^ = i^{^k)- Furthermore, the density of state, pk = p{^k) and the 
equilibrium populations p^'^ = p'^^{ek) scale as 1/A^. One thus can neglect a number of terms in 
the calculations. 

In the actual calculations, eq.(|4]) is used for the field-dependence of the tran sition rates. With 



this, one has for the relevant part of Lj^i{t2,ti) = Lj^i {t2 — h) according to eq.(A.5): 

L^^i\t) = e-^^'ni {pkXl - 4zXf) + Oil/N) 

where Xm = 'jMk — pMi and X/' = J2k PkX^i- The system is assumed to be in thermal equilibrium 
in the beginning and the field is assumed to be of the form H{t) = Hq cos (cot), cf. eq.([9]). The 
expressions given in Appendix A are used to compute the frequency-dependent response functions 



according to eq.(lO). For the variables Mk = M{ek), the choice discussed in the text is used, cf. 



eqns.(21) and (22). Furthermore, one can utilize detailed balance in the form 

Pk{fi) = HkPT with {K) = '^Kkpl'^ (B.l) 



k 



For the linear response one finds according to eq.(A.6) 



Xiicu) = {i + p)PY.pT(MI)^^^ (B.2) 

k Kk-lU 

which is just the Fourier transform of the time-derivative of the correlation function C^f (t). For 
the third-order response, one finds for n = 1, 3: 

Xt\^) = \P\1 + P) {xfe^(c^) + xgA(c^) + tU^) + Xg(^)} (B.3) 



where the individual terms are given by: 

xfeH^) = 3/i^EPfc(^')'4^i(^) -/^'EP^A(M5(Mf)5S(u;) (B.4) 

k k,l 

k,l k,l,m 



with 



-,(1) /■^^w ,_ _ / ^\'jl'^kf^l f^m ~^ ^ [oKkf^m ^l^k^l ^l^ll^r. 



^KSlil^{uj)) = Krr,Ki{K)- 






^ /c{3) / NN / .Kk^ll^m-^'^i'^l^k + ^l^l + Ql^m) 



-r /o(3) / NX / \ K,kK.l + 2KkK,jn + SKiKm - Q(^'^ 



This term corresponds to the first hne in Figjl} The second-order terms are: 

xgA(c-) = il-p) lpT.Pk{M'k) {^{MD - W)) S^aU^) (B.6) 



k,l 

where the corresponding spectral functions are: 



lT.Pm{Ml)UMl)-{M^))s'^ 






W 



Rpr c^W (, ^\ - Ik\ ^ 3/tfc^' + 2u^{AKk - t^i) 



Re(4;LM) = '^/('^)^- 



KfcKi — 6a;^ 



T /'c(3) / \\ / \1 2fi:fc + 3fi:i 



Additionally, I defined 



(M2)=Ep^(^') (B.8) 

k 



The other second-order term is: 



xgi^M = -3/^'Ep'c(M,2)^4;LM -7'(M2)Ep'c(^54;LH (B.9) 

k k 

-27/.Ep^P/(^5(M')4%M 
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with 






(3) ^, ,^,^ _ ,^ /,^\1 KkKi-3uj^ 

2{Kf + UJ^){Kl + 9uj^) 



R-el^gLM) = kz(k)^^_2 , 2^^ 2 






2 («:2 + a;2)(^2 + 9^2) 



Finally, the term corresponding to the third-order perturbation, i. e. the last line in FigjTj is 
given by: 

k k 

where 

Re(4^^M) = (/.)l ,,^^,, ; MS^^\co)) = {K)ly-^^ (B.12) 



2(4 + c^2) ' v-fe v-yy \ / 2 ^^ 



w^ 



Re(5f^M) = («:)i,-^^\— ; lm{SP{u)) = {^)- ^ 



Here, all expressions are given in a discrete notation. If one changes to a continuous description, 
one has to replace all sums by the appropriate integrals over the trap energies e^. 
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